'''
buildMesh.py - Paul Talbot, Oregon State University, NE 553, Spring 2011

Creates a mesh subdivided into regions, which are further subdivided into cells.
Regions are a group of cell with identical material properties.  Also creates
material properties which are subdivided into energy groups.  Each cell, region,
and mesh contains functions pertinent to a solution of the steady-
state neutron transport equation.  Used in conjunction with phiSolver.py,
runProblem.py, a quadrature file (extension .qdr), and an input file (extension .inp).
'''
import phiSolver as ps

class Mesh: #materials and subdivisions of problem
  def __init__(self,case=0,tolPhi=0.0001,quad=0):
    #attributes
    self.case='' #case being run
    self.ablys=[] #assemblies, groups of pins
    self.pins=[] #fuel pins, groups of regions
    self.rgns=[] #regions within mesh
    self.mtrls=[] #materials used in mesh
    self.mus=[] #angles of quadrature set
    self.wts=[] #weights that correspond to angles
    self.Fdx=0 #sum(Fission source * dx)
    self.Fco=[0 for i in range(5)] #coefficients for F
    self.K=1 #stores K-value
    self.tolPhi=tolPhi #phi convergence tolerance, used at cell level
    
    self.setQuadrature(quad) #establish quadrature

    self.Lcell=Cell(self,'dud','dud',0,0.1,0.0) #ghost cells used on boundaries
    self.Rcell=Cell(self,'dud','dud',0,0.1,0.0)
    self.Lrefl=0.0 #fraction reflecting
    self.Rrefl=0.0

    self.setCase(case) #set case to run
    self.readInput() #read input file, create mesh

    for rgn in self.rgns:
      rgn.setMaterial() #attach materials to regions

    self.setNeighbors()
    self.getFdx()

    for name in ['Uranium','MOX']:
      self.ablys.append(Assembly(self,name))

  def setCase(self,caseName):
    try:
      if caseName==0:
        caseName=raw_input('Enter case name: ') #accept user input
        #extension can be included or not
      if caseName.strip()[-4:]=='.inp':caseName=caseName[:-4]
      test=file(caseName+'.inp','r') #test to see if file can be opened
      test.close() #close test file so it can be opened later
      self.case=caseName #store case name as mesh attribute
      self.outFile=file(self.case+'.out','w')
    except IOError:
      print caseName+'.inp not found in current directory.\nPress Ctrl-D to exit.\n'
      self.setCase() #iterate back and try again

  def setQuadrature(self,quad):
    try:
      if quad==0: #no input quadrature given
        quad=raw_input('Enter quadrature input file: ') #get user unput
      quad=quad.strip() #remove extra spaces
      if quad[-4:]=='.qdr':quad=quad[:-4] #flexibility
      quadFile=file(quad+'.qdr','r') #try to open file
      #read in quadrature data
      flag=0 #1 angles, 2 weights
      for line in quadFile:
        if line[0]=='#' or line=='\n':
          continue
        elif line[:6]=='angles':
          flag=1
          continue
        elif line[:7]=='weights':
          flag=2
          continue
        elif flag==1:
          self.mus.append(float(line.strip()))
        elif flag==2:
          self.wts.append(float(line.strip()))
      if len(self.mus)!=len(self.wts):
        print 'Angles and weights are not evenly matched!  Exiting...'
        import sys
        sys.exit()
    except IOError:
      print quad+'.qdr not found in current directory.\nPress Ctrl-D to exit.\n'
      self.setQuadrature()

  def readInput(self):
    inFile=file(self.case+'.inp','r')
    flag=0 #gives the type of information being recorded, as follows:
      #0 None; 1 Region; 2 Materials; 3 Boundary Condition;
    for line in inFile:
      if line[0]=='#' or line=='\n':continue #skip empty or comment lines
      #if title line, change flag value
      elif line[:7]=='regions':flag=1
      elif line[:9]=='materials':flag=2
      elif line[:8]=='boundary':flag=3
      #if not a title line, store information
      elif flag==1: #region
        self.addRegion(line.split(';'))
      elif flag==2:
        if line[:8]=='material': #new material
          self.addMaterial(line[9:].strip()) 
        else:
          self.mtrls[-1].addGroup(line.split(';')) #more groups in material
      elif flag==3:
        if line.split(';')[0].strip()=='reflect':
          self.Lrefl=float(line.split(';')[1]) #boundary conditions
          self.Rrefl=float(line.split(';')[2])
        else:
          self.Lcell.phi=2*[float(line.split(';')[0])] #incident flux
          self.Rcell.phi=2*[float(line.split(';')[1])]

  def setNeighbors(self): #define neighbor cells, regions
    self.Lcell.setPos()
    for region in self.rgns:
      region.setNeighbors()
    self.Rcell.setPos()

  def addRegion(self,data): #add new region to mesh
    data[0]=data[0].strip() #name
    data[1]=float(data[1]) #width
    data[2]=int(data[2]) #number of cells
    data[3]=float(data[3]) #initial flux
    data[4]=int(data[4]) #fuel pin
    self.rgns.append(Region(self,len(self.rgns),data[0],data[1],data[2],data[3],data[4]))

  def addMaterial(self,name): #add material to mesh
    self.mtrls.append(Material(self,name))

  def printInfo(self): #useful for debugging
    print 'Mesh Attribute Output'
#    print '============='
#    print '| Materials |'
#    print '============='
#    for mtrl in self.mtrls:
#      mtrl.printInfo()
#    print '\n'
    print '=============='
    print '| Assemblies |'
    print '=============='
    for ably in self.ablys:
      ably.printInfo() #TODO write this function
    print '============='
    print '| Fuel Pins |'
    print '============='
    for pin in self.pins:
      pin.printInfo()
    print '\n'
#    print '==========='
#    print '| Regions |'
#    print '==========='
#    for rgn in self.rgns:
#      rgn.printInfo()

  def getFdx(self): #returns flux distribution function
    #get total F*dx
    Fdx=0
    for rgn in self.rgns: #update F*dx for all regions
      rgn.getFdx()
      Fdx+=rgn.Fdx
    self.Fdx=Fdx #store new value of F
    for rgn in self.rgns:
      for cell in rgn.cells:
        cell.F=cell.F/self.Fdx
    self.Fco[0]=self.Fdx

  def getK(self):
    self.K=self.Fdx*self.K

  def newGen(self,g):
    for rgn in self.rgns:
      rgn.newGen(g)

  def updateBoundaries(self):
    for g in [0,1]: #fast and thermal
      for n in range(len(self.mus)):
        if self.mus[n]>0:
          self.Lcell.g[g].psip[n]=self.Lrefl*self.rgns[0].cells[0].g[g].psim[n+5]
        else:
          self.Rcell.g[g].psim[n]=self.Rrefl*self.rgns[-1].cells[-1].g[g].psip[n-5]

  def doTransport(self,g,n):
    #do psi sweeps
    if self.mus[n]>0:
      for rgn in self.rgns: #sweep left to right
        rgn.doTransport(g,n)
    else:
      rvRgn=self.rgns[:]
      rvRgn.reverse() #sweep right to left
      for rgn in rvRgn:
        rgn.doTransport(g,n)
  
  def getPhi(self,g):
    check=True #all regions converged by default
    for rgn in self.rgns:
      rgnConv=rgn.getPhi(g)
      check=check and rgnConv #false if region not converged
    return check

  def writePhis(self,thermOut,fastOut):
    for rgn in self.rgns:
      for cell in rgn.cells:
        thermOut.writelines(str(cell.g[1].phi[0])+';')
        fastOut.writelines(str(cell.g[0].phi[0])+';')
    thermOut.writelines('\n')
    fastOut.writelines('\n')

class Assembly:
  class AblyGroup: #stores energy-dependent assembly data
    def __init__(self,ably,g):
      self.ably=ably #assembly which the group belongs to
      self.g=g #1 fast 2 thermal
      self.phi=0 #assembly-averaged phi
      self.phiCo=[0 for i in range(5)] #phi expansion coefficients
      self.sigt=0 #homogenized total xsec
      self.sigs_in=0
      self.sigs_out=0
      self.siga=0
      self.sigf=0
      self.sigr=0 #removal cross section, sig_t - sigs_in
      self.fL=0 #left discontinuity proportionality
      self.fR=0 #right
      self.D=0 #diffusion coefficient, 1/(3sigt)
    def avPhi(self):
      totPhi=0 #actually assembly total phi*dx
      totdx=0 #total assembly width
      for pin in self.ably.pins:
        totPhi+=pin.width*pin.groups[self.g-1].phi
        totdx+=pin.width
      self.ably.width=totdx
      self.phi=totPhi/totdx #FIXME needed or not?
      self.fL=self.ably.pins[0].rgns[0].cells[0].g[self.g-1].phi[0]/self.phi
      self.fR=self.ably.pins[-1].rgns[-1].cells[-1].g[self.g-1].phi[0]/self.phi
    def homogenize(self):
      t=0
      a=0
      f=0
      s_in=0
      s_out=0
      self.ably.nCells=0 #prevents pileup on multiple calls
      for pin in self.ably.pins:
        for rgn in pin.rgns:
          for cell in rgn.cells:
            self.ably.nCells+=1
            t+=cell.g[self.g-1].phi[0]*cell.mtrl.groups[self.g-1].sigt
            a+=cell.g[self.g-1].phi[0]*cell.mtrl.groups[self.g-1].siga
            f+=cell.g[self.g-1].phi[0]*cell.mtrl.groups[self.g-1].sigf
            s_in+=cell.g[self.g-1].phi[0]*cell.mtrl.groups[self.g-1].sigs_in
            s_out+=cell.g[self.g-1].phi[0]*cell.mtrl.groups[self.g-1].sigs_out
      self.sigt=t/(self.ably.nCells*self.phi)
      self.siga=a/(self.ably.nCells*self.phi)
      self.sigf=f/(self.ably.nCells*self.phi)
      self.sigs_in=s_in/(self.ably.nCells*self.phi)
      self.sigs_out=s_out/(self.ably.nCells*self.phi)
      self.sigr=self.sigt-self.sigs_in
      self.D=1.0/(3.0*self.sigt)
  #end class AblyGroup
  def __init__(self,mesh,name):
    self.mesh=mesh
    self.pins=[] #pin cells belonging to this group
    self.name=name
    self.mtrls=[] #list of materials in assembly
    for mtrl in self.mesh.mtrls:
      if mtrl.name[:3]==self.name[:3]:self.mtrls.append(mtrl)
    self.groups=[self.AblyGroup(self,1),self.AblyGroup(self,2)] #assembly group data
    self.nCells=0 #total number of cells in assembly
    self.setPins()
  def setPins(self):
    self.nCells=0 #prevents pileup over multiple calls
    self.nPins=0
    for pin in self.mesh.pins:
      for mtrl in self.mtrls:
        if mtrl in pin.mtrls:
          self.pins.append(pin)
          self.nCells+=pin.nCells
          self.nPins+=1
  def printInfo(self):
    print 'Material:',self.mtrl.name
    print 'No. Pins:',len(self.pins)
    print 'No. Groups:',len(self.groups)

class FuelPin: #subclass of mesh, bigger than region
  class PinGroup: #stores energy-dependent pin cell averages
    def __init__(self,pin,num):
      self.pin=pin #fuel pin to which this group belongs
      self.num=num #1 fast; 2 thermal
      self.phi=0 #pin average flux
  #end PinGroup class
  def __init__(self,mesh,num):
    self.mesh=mesh #mesh on which pin is found
    self.num=num #number of pin cell
    self.nCells=0 #number of cells in pin
    self.width=0 #width of pin cell
    self.groups=[self.PinGroup(self,1),self.PinGroup(self,2)] #stores fuel pin data by energy
    self.mtrls=[] #list of materials in pin
    self.rgns=[] #stores regions belonging to this group
  def getPhi(self,group):
    total=0 #total phi
    totdx=0 #total dx of cells
    nCells=0 #number of cells in pin
    for rgn in self.rgns:
      for cell in rgn.cells:
        nCells+=1
        total+=cell.g[group-1].phi[0]*cell.dx
        totdx+=cell.dx
    self.nCells=nCells
    self.width=totdx
    self.groups[group-1].phi=total/totdx #FIXME #weighted average
  def printInfo(self):
    print 'Pin Number:',self.num
    print '  Regions:',[rgn.num for rgn in self.rgns]
    print '  Groups:',[group.num for group in self.groups]
    print '  Materials:',[mtrl.name for mtrl in self.mtrls]
    print '  Width:',self.width
    print '\n'

class Material: #kinds of material cell can be made of
  def __init__(self,mesh,name):
    self.mesh=mesh #attaches it to mesh
    self.name=name #name of material
    self.groups=[] #energy groups to store data for

  def addGroup(self,data):
    name=data[0].strip() #group name (fast, thermal)
    sigt=float(data[1]) #total xsec
    sigs_in=float(data[2]) #in-group scattering xsec
    sigs_out=float(data[3]) #out-of-group scattering xsec
    siga=float(data[4]) #absorption xsec
    sigf=float(data[5]) #fission xsec
    nu=float(data[6])  #average neutrons per fission
    chi=float(data[7]) #fission distribution in energy
    self.groups.append(Group(self,name,sigt,sigs_in,sigs_out,siga,sigf,nu,chi))

  def findGroup(self,name): #find correct group, called by neutron class
    for group in self.groups:
      if group.name==name:
        myGroup=group
        break
    return myGroup

  def printInfo(self): #useful for debugging
    print 'Material Name: '+self.name+'.'
    print 'Number of Groups:',len(self.groups)
    for group in self.groups:
      group.printInfo()

class Group: #stores data for material dependent on energy level
  def __init__(self,mtrl,name,sigt,sigs_in,sigs_out,siga,sigf,nu,chi):
    self.mtrl=mtrl #material whose data is being stored
    self.name=name #fast, thermal
    self.sigt=sigt #material properties, as per above
    self.sigs_in=sigs_in
    self.sigs_out=sigs_out
    self.siga=siga
    self.sigf=sigf
    self.nu=nu
    self.chi=chi

  def printInfo(self): #useful for debugging
    s=self #shorthand
    print '  Group Name:',s.name,'; sigt:',str(s.sigt)[:3],'; sigsi:',\
        str(s.sigs_in)[:3],'; sigso:',str(s.sigs_out),'; siga:',\
        str(self.siga)[:3],'; sigf',str(self.sigf)[:3],'; nu:',str(s.nu)[:3],\
        '; chi:',str(self.chi)[:3]

class Region: #regions are groups of cells having same material
  def __init__(self,mesh,number,mtrl,width,nCells,initFlux,fuelPin):
    #attributes
    self.mesh=mesh
    self.num=number #ordial, from mesh build
    self.pinNum=fuelPin #number of fuel pin
    self.pin=0 #will be replaced by actual fuel pin instance
    self.materialName=mtrl #name of material
    self.mtrl='' #will be replaced by actual material instance
    self.width=width #length of region
    self.nCells=nCells #number of cells in region
    self.celldx=self.width/self.nCells #length of each cell
    self.initFlux=initFlux #initial flux in each cell in region
    self.Fdx=0 #total F*dx of all cells
    self.cells=[] #vector of cell objects

    #set fuel pin
    found=False
    for pin in self.mesh.pins:
      if self.pinNum==pin.num: #this region belongs to this pin
        self.pin=pin
        self.pin.rgns.append(self)
        found=True
        break
    if found==False:
      self.mesh.pins.append(FuelPin(self.mesh,len(self.mesh.pins)+1))
      self.pin=self.mesh.pins[-1]
      self.pin.rgns.append(self)

    #create cells
    for i in range(self.nCells):
      self.cells.append(Cell(self.mesh,self,self.mtrl,\
          len(self.cells),self.celldx,self.initFlux))

  def setMaterial(self): #attach actual material objects to region
    found=False
    for mtrl in self.mesh.mtrls:
      if mtrl.name==self.materialName:
        self.mtrl=mtrl
        found=True
        break
    if found==False:print 'Material not found for region:',self.materialName+'.'
    for cell in self.cells:
      cell.setMaterial()
    if self.mtrl not in self.pin.mtrls:
      self.pin.mtrls.append(self.mtrl)

  def addCell(self):
    self.cell.append(Cell(self.mesh,self,self.mtrl,self.len(cells),self.celldx,self.initFlux)) 

  def setNeighbors(self):
    #set left neighbor
    if self.num==0: self.lRgn='leftmost' #leftmost region
    else: self.lRgn=self.mesh.rgns[self.num-1] #all other regions
    #set right neighbor
    if self.num<len(self.mesh.rgns)-1: self.rRgn=self.mesh.rgns[self.num+1]
    else: self.rRgn='rightmost'
    #now go do it for the cells
    for cell in self.cells:
      cell.setNeighbors()

  def getFdx(self):
    Fdx=0
    for cell in self.cells:
      cell.getF()
      Fdx+=cell.F*cell.dx
    self.Fdx=Fdx

  def newGen(self,g):
    for cell in self.cells:
      cell.newGen(g)

  def doTransport(self,g,n):
    if self.mesh.mus[n]>0:
      for cell in self.cells:
        cell.doTransport(g,n) #one at a time, starting L going R
    else:
      rvCells=self.cells[:]
      rvCells.reverse() #get cells listed in reverse order
      for cell in rvCells:
        cell.doTransport(g,n) #one at a time, starting R going L

  def getPhi(self,g):
    check=True #all pass by default
    for cell in self.cells:
      cellConv=cell.getPhi(g) #sets phi and checks convergence
      check=check and cellConv #false if cell did not converge
    return check

  def printInfo(self): #useful for debugging
    print 'Region Number:',self.num
    print '  Material:',self.mtrl.name
    print '  Width:',self.width
    print '  L Cell Pos:',self.cells[0].pos
    print '  R Cell Pos:',self.cells[-1].pos
    print '  nCells:',self.nCells
    print '  celldx:',self.celldx
    print '  initFlx:',self.initFlux
    print '  cells:'
    for cell in self.cells:
      print 'L ',cell.lCell,'C ',cell,'R ',cell.rCell


class Cell: #smallest subdivided unit
  class CellGroup: #stores energy-dependent flux data for the cell
    def __init__(self,cell,num,initFlux):
      self.cell=cell #cell to which this group belongs
      self.num=num #1 fast; 2 thermal
      self.phi=[initFlux,initFlux,1] #phi current, old, older
      self.phix=[1,0,1] #phi slope current, old, older
      self.psip=[] #psi(i+1/2)
      self.psia=[] #psi(i) average
      self.psix=[] #psix(i) slope
      self.psim=[] #psi(i-1/2)
      self.Q=0 #fission source plus scattering source, value
      self.Qx=0 #slope of Q
      for psi in [self.psip,self.psim,self.psia,self.psix]:
        for n in self.cell.mesh.mus:
          psi.append(0) #forces psi to be as long as mus

  def __init__(self,mesh,rgn,mtrl,num,dx,initFlux):
    self.mesh=mesh #link to problem
    self.rgn=rgn #link to region
    self.mtrl=mtrl #link to material
    self.num=num #number within region
    self.dx=dx #width of cell
    self.pos=0.0 #length from L of cell to L of mesh, to be calculated
    self.lCell='' #cell to left, to be replaced
    self.rCell='' #cell to right, to be replaced
    self.F=1 #fission source term for cell, to be calculated
    self.Fx=0 #slope of fission source
    #create subclass to store group data
    self.g=[self.CellGroup(self,1,initFlux),self.CellGroup(self,2,initFlux)] 

  def setMaterial(self):
    self.mtrl=self.rgn.mtrl

  def setNeighbors(self):
    #define cell to left of this one
    if self.num>0: #not leftmost in region or mesh
      self.lCell=self.rgn.cells[self.num-1]
    elif self.rgn.lRgn=='leftmost': #leftmost cell in mesh
      self.lCell=self.mesh.Lcell
    else: #leftmost in region, not mesh
      self.lCell=self.rgn.lRgn.cells[-1]
    #define cell to right of this one
    if self.num<len(self.rgn.cells)-1: #not rightmost in region or mesh
      self.rCell=self.rgn.cells[self.num+1]
    elif self.rgn.rRgn=='rightmost': #rightmost in mesh
      self.rCell=self.mesh.Rcell
    else: #rightmost in region, not mesh
      self.rCell=self.rgn.rRgn.cells[0]
    self.setPos()

  def getF(self): #calculate F for cell
    self.F=self.mtrl.groups[0].nu*self.mtrl.groups[0].sigf*self.g[0].phi[0] +\
        self.mtrl.groups[1].nu*self.mtrl.groups[1].sigf*self.g[1].phi[0]
    self.Fx=self.mtrl.groups[0].nu*self.mtrl.groups[0].sigf*self.g[0].phix[0] +\
        self.mtrl.groups[1].nu*self.mtrl.groups[1].sigf*self.g[1].phix[0]

  def setPos(self): #calculate length from left side of mesh to left side of cell
    pos=0.0
    if self.rgn=='dud':
      if len(self.mesh.rgns)==0:
        self.pos=0.0
        return
      else:
        self.pos=float(self.mesh.rgns[-1].cells[-1].pos)+\
            float(self.mesh.rgns[-1].cells[-1].dx)
    else:
      i=0
      curRgn=self.mesh.rgns[i]
      while curRgn.num<self.rgn.num:
        pos+=curRgn.nCells*curRgn.celldx
        i+=1
        try:curRgn=self.mesh.rgns[i]
        except IndexError:
          #self.mesh.printInfo()
          #print curRgn.num,self.rgn.num
          print 'trying to get region',i
          import sys
          sys.exit()
      j=0
      curCell=self.rgn.cells[j]
      while curCell!=self:
        pos+=curCell.dx
        j+=1
        curCell=self.rgn.cells[j]
      self.pos=pos

  def newGen(self,g):
    #cycle storage variables
    self.g[g-1].phi=[0]+self.g[g-1].phi[:2]
    self.g[g-1].phix=[0]+self.g[g-1].phix[:2]

  def doTransport(self,g,n):
    #calculate constants
    tau=self.mtrl.groups[g-1].sigt*self.dx/abs(self.mesh.mus[n])
    import math as m
    A=m.exp(-tau)
    B=(1-m.exp(-tau))/self.mtrl.groups[g-1].sigt
    C=(1- 2/tau+ A*(1+2/tau))/self.mtrl.groups[g-1].sigt
    if g==1: #fast
      S=self.F/self.mesh.K #fission source
      Sx=self.Fx/self.mesh.K #fission slope
      T=self.mtrl.groups[g-1].sigs_in*self.g[g-1].phi[0] #in-scattering source
      Tx=self.mtrl.groups[g-1].sigs_in*self.g[g-1].phix[0] #in-scattering slope
    elif g==2: #thermal case
      S=self.mtrl.groups[0].sigs_out*self.g[0].phi[0] #downscatter source
      Sx=self.Fx/self.mesh.K #fission slope Sx
      T=self.mtrl.groups[g-1].sigs_in*self.g[g-1].phi[0] #in-scattering source
      Tx=self.g[g-1].phi[0]*self.mtrl.groups[g-1].sigs_in #in-scattering slope
    self.g[g-1].Q=0.5*(T+S) #Q calc
    self.g[g-1].Qx=0.5*(Tx+Sx) #Qx calc
    #pass to sweeping function
    if self.mesh.mus[n]>0:
      self.sweepRight(A,B,C,g,n) #from left to right
    elif self.mesh.mus[n]<0:
      self.sweepLeft(A,B,C,g,n) #from right to left

  def getPhi(self,g):
    #calculate new value
    self.g[g-1].phi[1:]=self.g[g-1].phi[:2] #store old values
    self.g[g-1].phix[1:]=self.g[g-1].phix[:2] #store old values
    self.g[g-1].phi[0]=0 #reset new value
    self.g[g-1].phix[0]=0 #reset new value
    for n,wt in enumerate(self.mesh.wts):
      self.g[g-1].phi[0]+=wt*self.g[g-1].psia[n]
      self.g[g-1].phix[0]+=wt*self.g[g-1].psix[n]
    #check tolerance on phi only
    check=ps.testTolerance(self.g[g-1].phi,self.mesh.tolPhi) #check phi tolerance
    return check

  def sweepRight(self,A,B,C,g,n):
    #calc psi(i+1/2)
    self.g[g-1].psip[n]=A*self.lCell.g[g-1].psip[n]+\
        B*self.g[g-1].Q + C*self.g[g-1].Qx
    #calc psi(i)
    self.g[g-1].psia[n]=(self.g[g-1].Q-self.mesh.mus[n]/self.dx*\
        (self.g[g-1].psip[n]-self.lCell.g[g-1].psip[n]))/self.mtrl.groups[g-1].sigt
    #calc psix(i)
    self.g[g-1].psix[n]=(self.g[g-1].Qx-3*self.mesh.mus[n]/self.dx*\
        (self.g[g-1].psip[n]+self.lCell.g[g-1].psip[n]-2*self.g[g-1].psia[n]))/self.dx

  def sweepLeft(self,A,B,C,g,n):
    self.g[g-1].psim[n]=A*self.rCell.g[g-1].psim[n]+\
        B*self.g[g-1].Q - C*self.g[g-1].Qx #calc psi(i+1/2)
    self.g[g-1].psia[n]=(self.g[g-1].Q-self.mesh.mus[n]/self.dx*\
        (self.rCell.g[g-1].psim[n]-self.g[g-1].psim[n]))/self.mtrl.groups[g-1].sigt #calc psi(i)
    self.g[g-1].psix[n]=(self.g[g-1].Qx-3*self.mesh.mus[n]/self.dx*\
        (self.rCell.g[g-1].psim[n]+self.g[g-1].psim[n]-2*self.g[g-1].psia[n]))/self.dx #psix(i)


# MAIN, used for debugging mesh builder
if __name__=='__main__':
  myMesh=Mesh()
  myMesh.printInfo()
